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expansions  of  the  equations  of  motion  as  the  frequency  of  the  forcing  function  tends  to  zero. 
We  characterize  such  interaction  as  being  generated  by  a  boundary  layer.  Such  a  layer 
justifies  the  discontinuity  that  certain  variables  must  undergo  across  this  interface  in  order 
to  numerically  integrate  the  Navier-Stokes  equations  with  boundary  conditions. 
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LOAD  NUMBERS,  SOLID  EARTH  TIDES, 
AND  LIQUID  CORE  DYNAMICS 


INTRODUCTION 

In  NRL  Report  8410  [1]  we  expressed  the  deformations  of  an  elastic  Earth  due  to  static  surface 
loads  as  infinite  series  of  spherical  harmonics.  The  coefficients  of  these  harmonics  were  called  the  load 
numbers  and  were  obtained  by  integrating  the  Navier-Stokes  equation  throughout  the  interior  of  a  lay¬ 
ered  Earth  with  boundary  conditions  that  had  to  be  satisfied  at  both  the  center  of  the  configuration  and 
the  loaded  upper  surface.  To  numerically  evaluate  these  infinite  series,  one  must:  (a)  be  able  to  ascer¬ 
tain  the  load  numbers  of  arbitrarily  high  order  without  having  to  solve  each  time  a  boundary  problem 
(this  is  tantamount  to  establishing  an  asymptotic  expansion  for  the  sequence  of  load  numbers)  and  (b) 
develop  practical  and  efficient  ways  for  numerically  summing  those  series  of  spherical  harmonics  once 
the  load  numbers  are  known. 

We  discuss  in  this  report  a  question  that  was  left  unanswered  in  our  previous  publication,  that  is: 
why  the  stress  function  was  assumed  to  be  discontinuous  at  the  core-mantle  interface  when  numerically 
evaluating  the  load  numbers  as  a  boundary  value  problem. 

This  report  presents  in  a  systematic  and  rigorous  way  procedures  suitable  for  the  numerical 
evaluation  of  the ^Earth’s  spheroidal  deformations,  which  most  commonly  are  called  the  Earth’s  tides. 
For  this  purpose,  we:  (a)  study  the  Boussinesq  theory  for  the  elastic  deformations  of  a  flat  plate  and 
show  how  by  applying  the  Boussinesq  solution  to  the  tangent  plane  at  a  loading  point  of  the  spheroidal 
Earth  we  can  establish  the  asymptotic  behavior  of  the  load  numbers  [2];  (b)  elaborate  on  certain 
numerical  procedures  which  seem  to  be  best  suited  for  the  summation  of  series  once  the  asymptotic 
behavior  of  their  coefficients  has  been  ascertained;  (c)  provide  in  tabular  form  the  closed  form  expres¬ 
sions  of  those  infinite  series  of  spherical  harmonics  which  are  essential  for  the  numerical  evaluation  of 
the  Earth’s  tides;  and  (d)  examine  the  dynamics  of  a  liquid  core  for  a  nonrotating  Earth  to  mathemati¬ 
cally  prove  the  existence  of  a  boundary  layer  at  the  top  of  the  liquid  core.  This  layer  provides  a 
justification  for  the  discontinuous  behavior  of  some  of  the  integration  variables  at  the  liquid  core/solid 
mantle  interface,  and  furnishes  the  proper  number  of  free  parameters  for  satisfying  the  boundary  condi¬ 
tions  at  the  loaded  surface. 

/K 

A  more  detailed  discussion  of  these  ideas  can  be  found  in  a  forthcoming  publication  by  Lanzano 
[31.  In  reaching  our  main  results,  we  want  to  acknowledge  the  fact  that  we  have  greatly  benefited  from 
two  older  publications  by  Farrell  [4]  and  Pekeris  and  Accad  [51. 

ASYMPTOTIC  VALUES  OF  THE  LOAD  NUMBERS 

The  strategy  we  follow  in  order  to  ascertain  the  asymptotic  values  of  the  load  numbers  can  be 
summarized  in  the  following  steps: 

(I)  One  should  compare  the  infinite  series  expansions  which  represent  the  spheroidal  defor- 
.  mations  and  which  contain  the  load  numbers  as  the  coefficients  of  the  various  harmonics 
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with  appropriate  closed  form  solutions  for  the  same  type  of  deformations  which  are  valid 
within  a  limited  region  of  the  spheroid.  This  is  the  case  in  the  neighborhood  of  a  loading 
point  where  the  Boussinesq  theory  of  the  flat  plate  applies  and  gives  rise  to  a  closed  form 
solution. 

(2)  One  should  refer  both  solutions  in  the  region  where  they  overlap  to  a  common  reference 
frame.  It  is  convenient  to  make  use  of  the  cylindrical  coordinate  system  with  the  loading 
line  as  its  symmetry  axis  and  the  tangent  plane  to  the  spheroid  at  the  loading  point  as  the 
auxiliary  plane. 

(3)  Because  the  limit  of  the  Legendre  polynomials  Pn  for  large  values  of  n  is  the  Bessel  func¬ 
tion  Jo  (of  order  zero)  [6],  and  similarly  the  same  limit  for  the  derivative  dPj d9  is  pro¬ 
portional  to  J{,  one  should  solve  the  Boussinesq  problem  by  means  of  the  Hankel 
transforms  of  order  zero  and  order  one  because  these  functional  transforms  use  those 
Bessel  functions  as  their  kernels. 

(4)  As  a  consequence,  it  follows  that  the  asymptotic  values  of  the  three  load  numbers  must 
be  proportional  to  the  Hankel  transforms  of  the  two  components  of  the  displacement  and 
of  the  perturbed  potential  in  the  Boussinesq  problem. 

We  shall  next  elaborate  statements  1  through  4.  To  begin  with,  it  is  clear  that  for  small  angular 
separations  from  the  loading  line,  the  displacements  and  the  perturbed  potential  of  a  layered  spheroid 
must  agree  with  the  corresponding  quantities  referring  to  the  tangent  plane  at  the  loading  point.  In  the 
neighborhood  of  a  loading  point,  and  for  the  values  of  density  and  elastic  parameters  commonly 
accepted  in  the  most  recent  Earth  models,  the  elastic  forces  dominate  the  gravitational  forces.  We  can, 
therefore,  not  only  assume  that  the  density  and  Lamd  parameters  remain  constant  in  the  neighborhood 
of  a  loading  point,  but  also  neglect  the  coupling  between  the  gravitational  and  elastic  forces;  this  situa¬ 
tion  enables  us  to  study  and  separately  solve  one  elastic  problem  and  one  gravitational  problem. 

In  the  absence  of  the  gravity  field  g0  of  the  reference  state  and  of  the  perturbed  gravity  gj  due  to 
deformation,  assuming  also:  (1)  elastic  equilibrium,  and  (2)  constancy  of  the  Lam6  parameters,  the 
Navier-Stokes  equation,  which  was  developed  in  our  previous  work,  Eqs.  (14)  and  (15)  of  NRL  Report 
8410  [1],  will  simply  reduce  to 

(X  +  2/4.) V (V  •  V)  -  mV  x  V  x  V,  (1) 

where  u  is  here  the  displacement  vector.  The  above  equation  should  be  solved  for  the  material  half¬ 
space  simulating  the  tangent  plane,  under  the  assumption  that  there  is  only  one  normal  load  concentrat¬ 
ed  at  a  given  location  in  the  uppermost  surface.  By  the  same  token,  the  perturbations  in  the  gravita¬ 
tional  potential  can  be  represented  by  the  Poisson  equation 

V2<*>  -  -4irGp„(V  ■  17).  (2) 

with  constant  density  p0,  along  the  material  half-space,  and  by  the  Laplace  equation 

V2*  -  0, 

in  empty  space. 

Consider  a  spheroidal  Earth  bounded  by  an  equipotential  surface _r  —  a  with  a  vertical  load  applied 
at  one  of  its  points  Q.  We  introduce  the  spherical  reference  frame  (OQ,  r,  9)  with  origin  Oat  the  cen¬ 
troid  of  the  Earth  and  polar  axis  along  the  line  OQ,  we  shall  measure  the^olatitude  9  from  this  axis  and 
assume  that  the  configuration  be  symmetrical  with  respect  to  the  line  OQ.  This  line  of  application  of 
the  load  can  be  taken  as  the  r-axis  of  a  cylindrical  reference  with  origin  at  the  loading  point  Q,  and  the 
tangent  plane  to  the  spheroid  at  Q  can  be  taken  as  the  auxiliary  plane  upon  which  the  other  two  vari¬ 
ables  (rand  *)  will  be  measured.  This  cylindrical  system  shall  be  denoted  by  (Q,  r,  <li).  It  is  clear  that 
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if  we  take  a  point  T  in  the  neighborhood  of  Q,  the  arc  (QT)  along  the  equipotential,  which  is 
represented  by  a9 ,  must  be  the  radial  distance  in  the  cylindrical  reference  system;  also  that  the  unit 
vectors  ez  and  er  of  the  cylindrical  frame  must  be  the  limiting  positions  of  the  unit  vectors  e,  and  e», 
respectively,  of  the  associated  spherical  frame. 

We  shall  solve  both  Eqs.  (1)  and  (2)  where  the  displacement  vector  V  will  have  the  following 
representation  in  cylindrical  coordinates: 

V  -  u(,2,r)ez  +  v(z,r)er.  (3) 

The  spheroidal  deformations  of  a  layered  spheroid  along  an  equipotential  surface  r  —  a  due  to  an 
applied  vertical  load  of  mass  mo  acting  at  a  point  Q,  and  when  expressed  in  spherical  coordinates  in 
terms  of  the  load  numbers  h\  /',  can  be  written  as  follows: 

Via, 9)  =  T  [e,  h'„{a)P„(cos  9)  +  ee  l‘„{a)dP„{cos  9)/d9],  (4) 

m  n-0 

Here  m  is  the  mass  of  the  spheroid  within  the  equipotential  surface  r  —  a. 

Similarly,  the  potential  due  to  the  deformation  of  the  spheroid  can  be  expressed  in  terms  of  the 
third  load  number  k\  also  in  spherical  coordinates,  as  follows: 

<t>(a,9)  - -  y0 (a)  T  k„(a)Pn( cos  9).  (5) 

m  n-0 

Here  yo(a)  is  the  gravitational  attraction  since  we  are  neglecting  the  rotational  effects. 


It  is  now  appropriate  to  recall  that  the  limit  of  the  Legendre  polynomials  for  large  values  of  the 
parameter  n  can  be  established  in  terms  of  the  Bessel  functions.  More  specifically,  see  Ref.  6,  p.  1 55 


lim  P n [cos  (9/n)]  -  J0(9).  (6) 

IT— oo 


By  differentiation,  and  if  use  is  made  of  elementary  properties  of  Bessel  functions, 

-J,(9). 


1  dP„  [cos  (0/n)] 
-I™  n  dWn ) 


(7) 


Because  of  these  properties  it  is  convenient  to  solve  Eqs.  (1)  and  (2)  by  means  of  the  Hankel 
transforms. 


The  Hankel  transform  of  order  n  of  a  function  f(z,  r)  of  two  variables  with  respect  to  one  of  its 
variables  is  represented  by  the  corresponding  capital  letter  and  is  defined  as 

F(z,  *)  -  fQ  fiz,  r)J„{(r)rdr,  (8) 

where  J„  is  the  Bessel  function  of  order  n.  We  can  verify  that  the  inversion  of  Eq.  (8)  will  provide 

/fc  r)-f~F(z.  ()JMr)(d(.  (9) 

Note  that  the  product  (r  of  these  two  variables  must  be  a  nondimensional  quantity.  From  Eqs. 
(8)  and  (9),  it  is  clear  that  the  dimensions  of  a  set  of  Hankel  transforms  (/;/")  will  differ  by  the  square 
of  the  dimension  of  the  f  or  r- variable. 

We  now  revert  to  Eq.  (1)  and  write  the  components  of  the  stress  as 

ra  -  A(e„  +  <M  +  fB)  +  2pfB  -  (A  +  +  ~ 

He 
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Here  u ,  v  are  the  components  of  the  displacement  vector  17 given  by  Eq.  (3);  the  notation  we  use  is  the 
one  by  Sokolnikoff  [7]. 

Consider  the  system  formed  by  Eqs.  (1)  and  (10);  and  apply  the  Hankel  transform  of  order  zero 
to  the  u  and  za  variables,  the  H  transform  of  order  one  to  the  v  and  ra  variables.  Denote  the  H 
transforms  of  these  four  variables  by  the  corresponding  capital  letters.  Then  assume  that  both  u  and  v 
and  their  partial  derivatives  vanish  at  infinity  with  Mr.  Consider  the  system  of  equations  in  the 
material  region  of  the  linear  space,  which  we  assume  to  be  the  z  <  0  half  space. 

We  reach  the  linear  differential  system 

T-  AS, 

where  the  primes  denote  derivatives  with  respect  to  the  ^variable  and  where  we  have  set 

3?  s  Column  (U,  V,  Ta,  Ta). 

The  matrix  A  is 

0  — Xf/<r  1/cr  0 

*  o  oi/m 
0  0  0  -£ 

0  A/j.r)^/cr  Xf/o-  0 

with  constant  values  for  <r  -  X  +  2  m  and  tj  -  X  +  ft. 

The  general  solution  of  Eq.  (11)  is  obtainable  by  elementary  procedures  in  terms  of  the  eigen¬ 
values  p  of  the  A  matrix,  i.e.,  in  terms  of  the  roots  of  the  characteristic  equation 

Det  04  -  pi)  -  0,  (14) 

where  I  here  stands  for  the  (4  x  4)  identity  matrix.  The  four  roots  of  Eq.  (14)  are  0i  -  Pi  -  f  and 
03-04--*. 


(11) 

(12)  \ 


The  eigenvectors  or  solutions  to  Eq.  (11)  can  be  expressed  according  to  the  infinite  relationship 


3?  -  exp  (±fz)  exp  [(A  =F  (I)z]y  -  exp  (±fz)l?  +  z(A =F  f  I)y  +  y(/4  =F  f/)2?  +  . 

..].  (15) 

One  can  obtain  two  independent  solutions. 

"V 

exp  (±{ z)7x, 

by  solving  the  linear  system. 

(16) 

(A  T  £/)?,  -  0. 

(17) 

Two  additional  independent  solutions  are 

l 

exp  (±{z)172  +  zC4  Tf/)y2], 

(18) 

i 

obtainable  by  choosing  solutions  of 

i 

i 

U  T  (I)1?!  -  0, 

(19) 

) 

which  are  not  common  to  Eq.  (17).  In  both  cases,  the  infinite  expansion  appearing  in  Eq.  (IS)  will  ter- 
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By  following  the  preceding  procedure,  we  reach  the  set  of  four  linearly  independent  solutions: 


1 

=fl 

±2 fi£ 

-2 #if 


exp  (±fz); 


(20) 


Tfz  +n/V 

fz  ±  o7t) 
"2 H(2z 
±2n$2z  +  2  fi( 


exp  (±fz). 


Let  us  now  take  into  account  the  boundary  conditions.  The  vanishing  of  the  four  original  variables  at 
z  —  — oo,  means  that  the  transformed  variables  U,  V,  Ta,  Tn  must  vanish  at  minus  infinity.  If  we 
choose  henceforth  f  >  0,  then  the  solutions  containing  exp  (-fz)  cannot  be  taken  into  consideration. 
To  combine  the  remaining  two  solutions,  we  suppose:  (a)  the  transversal  component  of  the  stress  van¬ 
ishes  at  z  —  0:  rn(0,  r)  —  0;  and  (2)  the  unit  normal  load  is  applied  on  a  spherical  cap  of  radius  r  — 
a,  i.e., 

ra(0,  r)  »>  —  1/w  a2  for  r  <  a, 
ra(0,  r)  -  0  for  r  >  a. 

It  follows  then  that 


and 


TJ0.&-0. 


Tb(0,() - l—  C  JQ{ri)r  dr 

it  a  *  u 


1  2  J\(£a) 

2l r  fa 


whose  limit  for  a  shrinking  a  is  —  —  l/2ir.  Using  the  previous  relations,  we  finally  have 


U(z.O 


1 


exp  (fz); 


K(z,f) 


1 


fz  +  exp  (fz). 


(21) 


We  next  consider  the  gravitational  problem  by  solving  the  Poisson  equation  within  the  material  tangent 
plane  (z  <  0)  and  the  Laplace  equation  in  empty  space  (z  >  0). 

We  introduce  the  Hankel  transform  of  order  zero  for  the  perturbed  potential  4>  and  denote  it  by 
d>.  The  Poisson  equation  then  becomes 


<Mz,f) 


2<jpo 

V 


exp  (fz). 


(22) 
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after  use  has  been  made  of  Eq.  (21)  to  express  the  components  U  and  V  of  the  displacement.  The 
characteristic  exponents  of  Eq.  (22)  are  ±£.  The  solution  of  the  Laplace  equation  valid  for  empty 
space  (z  >  0)  must  vanish  at  z  =  +°°  and  must  therefore  be  in  the  form 

<J>,(z,f)  —  A  exp  (—  |z).  (23) 


The  solution  of  Eq.  (22),  valid  within  the  elastic  medium,  can  be  seen  to  be: 


B  + 


Gpq  | 

vi  1 


exp  (fz). 


(24) 


However,  must  be  continuous  at  the  interface  z  —  0;  this  leads  to  A  —  B.  By  rewriting  the  Poisson 
equation  as  the  divergence  of  the  vector 


V0  +  4w  GpoV, 

we  realize  that  the  normal  component  of  such  vector  must  be  continuous  at  z  -  0.  In  terms  of  the  H 
transforms,  we  reach  the  condition: 


ld4>e/dz lj_o”  [d<J\/3z]2_o  +  AnGpoUjiO.g): 


(25) 


this  is  so  because  of  the  vanishing  of  Ut(z,£).  Use  of  Eq.  (25)  is  instrumental  in  determining  the 
value  of  A.  We  get 


<Mz.f) 


Gp(\ 

- \  exp  (-fz)  for  z  >  0,  and 


*/Uf) 


Cpo 

2pi2 


1  + 


exp  (f  z)  for  z  <  0. 


(26) 


We  are  now  in  the  position  of  comparing  two  solutions  which  are  both  valid  in  the  neighborhood  of  a 
loading  point  and  of  determining  the  limit  of  the  three  load  numbers  as  the  order  of  the  harmonic 
tends  to  infinity. 


First,  by  expressing  the  solutions  given  by  Eqs.  (4)  and  (5)  in  cylindrical  coordinates  and  by  mak¬ 
ing  use  of  Eqs.  (6)  and  (7),  the  displacement  vector  and  the  perturbed  potential  for  large  values  of  n, 
may  be  written: 


— -  lezh'(a)Jo(n9)  -  e,nl'„(a)J\(n9))\ - -  ya(a) kn(a)J0(nO) .  (27) 

m  m 

We  next  represent  the  Boussinesq  solution  for  z  =  0  and  for  small  values  of  the  radial  distance  r  —  a9 
as  integrals  of  the  H  transforms  according  to  the  following  relationships: 


S 1(0,a9)  -  e:  f“  i/jo.-j] /<,(«•) -jr  dn  +  e,  J"  pjo.-jJ./, (»*)—-  dn\ 
*(O,a0)  -  J0”  $  0,^-j  Jo(,n0)-^  dn. 


(28) 


Note  that  the  integration  variable  f  has  been  relabeled  f  —  n/o,  a  legitimate  replacement  since  f  has 
the  dimension  of  an  inverse  length. 
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Each  integral  appearing  in  Eq.  (28)  can  be  visualized  as  an  infinite  summation;  their  terms  per¬ 
taining  to  large  values  of  n  must  coincide  with  the  corresponding  terms  appearing  in  Eq.  (27).  Since 
the  Boussinesq  solution  was  written  in  terms  of  a  unit  force  whereas  the  applied  load  here  is  may0(a), 
we  reach  the  following  relationships: 

hn(a)  “  P" 

— -  nl'„(a)  -  ~  -Tf  I'jo,— I m0y0(a)\  (29) 

— ~  y0(a)k'„(a)  -  — y  4>|o,— |m0yo(a). 
m  a1  \  a j 

By  using  Eqs.  (21)  for  U  and  V,  and  Eqs.  (26)  for  <!>,,  we  finally  arrive  at  the  asymptotic  representation, 


A* 

/* 

myo 

1 

A* 

Aira1^ 

~  3p0T|/2jx,p 

Here  a  —  A  +  2/a,  17  —  X  +  /x,  p  is  the  mean  density,  m  is  the  total  mass  of  the  spheroid,  and  A*,  /*,  k * 
are  the  limits  of  h„,  nl'„,  nk'„  as  n  approaches  infinity.  All  quantities  appearing  in  Eq.  (30)  must  be 
evaluated  at  the  outermost  surface,  r  —  a. 

For  the  1066 A  Earth  model,  these  limits  can  be  calculated  and  are 

A*  **  —  11.35767, 

/*  -  3.43096,  (31) 

A*  «  -4.70473. 

NUMERICAL  EVALUATION  OF  THE  GEOPHYSICAL  PARAMETERS 

Once  the  asymptotic  behavior  of  the  three  load  numbers  has  been  ascertained,  one  can  numeri¬ 
cally  evaluate  the  elements  of  geophysical  interest;  these  are:  (a)  the  elastic  deformations  of  the  sur¬ 
face;  (b)  the  variation  of  the  gravity  field  in  direction  and  intensity;  and  (c)  the  strain  tensor  induced 
by  the  deformations.  All  these  quantities  are  represented  by  infinite  series  of  Legendre  polynomials, 
the  angular  distance  9  being  measured  from  the  loading  line. 

Because  we  know  the  asymptotic  limits  A*,  A*,  /*  of  h‘„,  nk„ ,  we  can  take  advantage  of  the 
Kummer  method  for  the  summation  of  an  infinite  series  of  functions  18,9).  This  method  consists  of 
subtracting  from  the  given  series  an  appropriately  chosen  series  of  known  sum  whose  elements  are 
asymptotically  proportional  to  the  series  in  question.  The  result  can  then  be  expressed  as  the  sum  of 
two  infinite  series  of  Legendre  polynomials:  (a)  the  first  one  is  multiplied  by  the  known  asymptotic 
value  of  the  load  number  and  its  sum  can  be  represented  in  closed  form;  (b)  the  second  infinite  series 
has  coefficients  asymptotically  tending  to  zero,  so  that  it  can  ultimately  be  evaluated  as  a  finite  sum;  the 
number  N  of  its  terms  must  be  so  chosen  that  for  values  of  the  subscript  n  larger  than  N,  the  difference 
between  the  corresponding  load  number  and  its  asymptotic  value  can  be  considered  negligible  according 
to  the  degree  of  approximation  we  plan  to  achieve. 


Let  us  further  examine  these  ideas.  Consider  the  components  u,  v  of  the  displacement  vector 
V  (a, 6).  For  the  normal  component,  we  get 


uUO) i  kHa)P„  B  (a*  i  Pn  +  £  <*»’  -  h')P*  ■ 


(32) 


the  value  of  N  depending  on  the  accuracy  to  be  attained. 
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In  a  similar  manner,  the  tangential  component  v  of  the  displacement  can  be  evaluated  according 
to  the  scheme: 


v(a,0 ) 


am0 


/*  T  — 

n 


I 

n  d9 


I  «' 


/•)1^ 


(33) 


This  is  so  because  /*  is  the  limit  of  «/„'  as  n  goes  to  infinity.  Both  infinite  series,  P„  and 
(1/n)  ( dPJdQ ),  that  appear  in  Eqs.  (32)  and  (33)  can  be  summed  in  closed  form;  their  expressions  are 
given  in  Table  1. 


Table  1  —  Closed  Form  Sum  of  Certain  Harmonic 
Series  of  Geophysical  Interest 


(1)  V  P„( cos  0)  -  y  cosec  ~ 
n-o  2  2 


(2)  T  nP„  (cos  0)  =  -  -j  cosec  I 

-Tn  4  |  2 

,2 


,31 


cosec 


(4)  £ 

n—  1 
oo 

(5)  £ 

f»-I 

oo 

(6)  £ 

1 


I 

n  dO 

1  ^ 


cos  9  +  sin  |y 


cosec0 


n  d9 2 

—  P„(  cos  0) 
n 


1  -  sin3  ^ 1|  cosec20 


-In 


sin 

9_ 

...  0 

1  +  sin  — 

2 

2 

Next,  consider  the  variation  in  the  gravity  field.  The  potential  at  a  point  Q(a  +  u,  9)  of  the 
deformed  surface,  located  at  a  height  u  from  the  point  P(a,9)  on  the  equipotential  r  —  a,  can  be  writ¬ 
ten: 

Y(Q)  =  F0(P)  +  uVo(P)  +  K*(/>)  +  V'(P). 

Here,  primes  denote  derivatives  with  respect  to  the  normal  distance,  V0  is  the  potential  of  the  original 
configuration,  F*  is  the  potential  due  to  the  redistribution  of  mass,  and  V,  is  the  potential  of  the 
applied  load. 

A  gravimeter  located  at  Q  measures  the  following  acceleration: 

V(Q)  -  Vq  ( p )  +  uVq(P)  +  ( V’)'{P )  +  y;(P).  (34) 

We  now  expand  both  F*  and  V,  into  spherical  harmonics  and  note  that  Vq  (P)  —  -yo (a),  where 

y0(a)  -  J0  p0(r)r2dr.  (35) 

By  differentiation,  we  easily  establish  that 

yo (a)  "  —  —  y<>(a)  +  47rGpo(fl).  (36) 

a 
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If  we  neglect  the  contribution  due  to  the  density  at  P ,  Eq.  (34)  can  be  rewritten  as 

-y(0)  -  -y0  (P)  +  ~  uy0(P)  -  —  V'X(P)  +  -  K,(/>), 
a  a  a 

which  leads  to 

y'iP.Q)  -  y0(F)  “  y<<?)  -  —  y0(/*>  £  [2*; (a)  -  (a  +  1)*>)  +  «]!»,.  (37) 

Finally,  the  tilt  or  deviation  of  the  vertical  direction  can  be  expressed  as  the  angle  between  the  normal 
to  equipotential  surface,  V0  +  V‘  —  constant  and  the  geometrical  normal  to  the  deformed  surface.  For 
simplicity,  we  shall  limit  our  considerations  to  the  component  of  the  tilt  upon  the  meridian  of  the 
spheroid  which  we  denote  by  /*  (a,  0). 


The  normal  line  to  the  equipotential  surface  is  represented  by 

(1  +  O(dPjdO), 

whereas  the  normal  to  the  deformed  surface  depends  on 

(1  /a)(dr/dO), 

which  is  proportional  to 

hi(dPJde). 

We  can  then  write 


t^a.O)  -  —  f  (1  +  *„’(«)  -  A»)  (38) 

m  do 

The  numerical  evaluation  of  Eqs.  (37)  and  (38)  requires  the  sum  of  three  additional  infinite  series. 
nP„ ;  (1  /n)P„\  dPjdO.  Their  sums  are  also  available  from  Table  1. 


We  now  come  to  the  evaluation  of  the  components  of  the  strain  tensor.  We  use  the  fundamental 
Navier-Stokes  equation  which  was  expressed  as  Eq.  (39)  in  Ref.  1  to  write: 


«„(a,0)  »  u’fa.O)  —  £  U„(a)Pn 
0 


I 

0 


2A  (a)  ,r  /  \  i  /  ,  \(u) 

UAa)  +  n(n  +  1)  - 


r(a) 


a<r(a) 


W  + 


E„(a) 


r(o) 


P.  (cos  0). 


As  usual,  primes  denote  derivatives  with  respect  to  the  radial  distance  and  cr  —  A  +  2/u..  U,  V  and  E 
denote  the  expansions  into  spherical  harmonics  of  the  components  «,  v,  of  the  displacement  vector  and 
of  the  component  of  the  stress  tensor.  We  know  that  E„(a)  is  a  delta- function  with  its  peak  at  0  - 
0,  so  that  it  will  give  a  zero  contribution  anywhere  else.  If  we  take  these  facts  into  account  and 
represent  the  v  component  of  the  displacement  vector  in  terms  of  the  load  numbers,  we  can  rewrite  the 
previous  equation  as 


e"(a’9)  -  ~  u(a’0)  +  ?  i nin  +  l)'»'(fl)/>«(cos  9)- 


We  can  next  verify  that 

1 

- 


u(a,0)  +  —  v(a,l0) 


1  ,  ,  Wn  “  ,,  ,  d*P. 

-  S  %  <■“  ST- 


(39) 


(40) 


Also 


—  lu(a,  0)  +  v(a,0)  cot  0], 
a 
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The  other  components  ot  the  strain  tensor  can  be  shown  to  vanish  when  evaluated  at  the  surface  r  -  a. 
The  numerical  computation  of  these  infinite  series  which  represent  the  components  of  the  strain  tensor 
necessitates  the  knowledge  of  the  sum  of  the  series  (1  /  n){d2Pj  d02)\  t’  :s  is  also  given  in  Table  1. 

Table  1  provides  the  closed  form  sum  of  all  those  series  which  are  required  for  the  numerical 
evaluation  of  the  geophysical  parameters.  A  detailed  derivation  of  the  sum  of  these  series  appears  in 
Ref.  3;  for  some  partial  results  see  Ref.  10. 

From  a  computational  point  of  view,  few  remarks  are  in  order  before  bringing  this  section  to  a 

close. 


(1)  The  second  infinite  series  that  was  introduced  through  the  Kummer  transformation  turns  out  to 
be  slowly  convergent;  various  numerical  algorithms  must  be  used  in  order  to  accelerate  its  conver¬ 
gence;  in  this  regard,  the  Euler  transformation  was  found  to  be  a  useful  tool  [8,9]. 

(2)  For  large  values  of  the  order  n  of  the  harmonic  and  for  small  values  of  the  angle  0  measured  from 
the  loading  line,  one  can  achieve  computational  economy  by  approximating  P„(cos  0)  and  dPJdB 
by  means  of  J0  and  J\,  as  has  been  described  in  the  previous  section. 

(3)  The  functions  appearing  in  Table  1  are  not  finite  at  8  —  0  and  should  not  be  used  in  the  immedi¬ 
ate  neighborhood  of  the  loading  position.  We  have  remedied  to  this  situation  by  employing  the 
Boussinesq  approximation  extended  to  a  suitable  neighborhood  of  the  loading  point.  We  resume 
therefore  our  considerations  on  the  Boussinesq  solution  which  had  reached  the  representation 
given  by  Eq.  (21)  and  try  to  express  it  in  terms  of  cylindrical  coordinates  at  the  loading  point.  In 
other  words,  we  shall  invert  the  Hankel  transformation.  For  this  purpose,  we  use  the  following 
integral  appearing  in  Ref.  6: 

J0  exp(-<jz)  Jy(bz)dz  «•  (C/bYia2  +  b2)~'/2, 
which  is  valid  for  v  >  0,  a  0,  b  ^  0  and  where  C  —  ( a 2  +b2)i/2  -  a  >  0. 


Differentiation  of  the  above  expression  with  respect  to  the  parameter  a  yields 

J"fl  exp  (—  az)Jy(bz)z  dz  —  [pb(C/b)v~l  +  a (.  1  —  p)(C/b)*](.a 1  +  b2)~2/2. 

The  preceding  two  integral  relations  for  v  —  0,  I  are  all  that  is  required  to  obtain  the  Boussinesq 
approximation  in  terms  of  cylindrical  coordinates  (z,  r): 


u(z,  r ) 


m0y0(a) 

4ir/iR 


(42) 


v(z,  r) 


moyp  (a)  .  j_  t )r2z 

4 irrjr  R  M R3 1 


Here  we  have  set 


R2  —  r2  +  z2;  rj  —  X  +  n  and  <r  —  X  +  2 fi. 


LIQUID  CORE/SOLID  MANTLE  DYNAMICS 

To  justify  certain  assumptions  that  were  made  in  carrying  out  the  numerical  integration  of  the 
Navier-Stokes  equations,  it  is  appropriate  to  examine  the  physical  conditions  that  might  be  expected  at 
the  interface  between  liquid  core  and  solid  mantle.  We  shall  therefore  revert  to  the  linearized  version 
of  the  Navier-Stokes  equations  which  was  obtained  in  our  previous  work  [1]  and  establish  some  funda¬ 
mental  properties  of  its  solutions. 
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For  the  scope  of  a  first-order  approximation,  we  can  safely  ignore  the  influence  of  the  Earth’s 
rotation;  however,  we  want  to  study  those  solutions  of  the  Navier-Stokes  equations  that  represent  har¬ 
monic  oscillations  depending  on  the  time  factor  exp  ( i<rt ).  We  believe  that  a  deeper  understanding  can 
be  obtained  of  the  conditions  existing  at  the  subject  interface  if  we  consider  the  Earth’s  permanent 
deformations  as  the  limit  of  spheroidal  oscillations  of  the  harmonic  type  when  their  frequency  <r  tends 
to  zero. 

We  use  here  the  same  notation  adopted  in  the  previous  report;  in  particular,  recall  that  the  three 
sets  of  functions  U. ,  K  and  R  represent  the  coefficients  of  the  spherical  harmonic  developments  for  the 
components  u,  v  of  the  displacement  vector,  and  for  the  potential  due  to  deformation. 

Within  the  liquid  core  where  p  “  0  and  if  we  ignore  the  effects  due  to  rotation,  the  equations 
governing  the  spheroidal  oscillations  of  constant  frequency  cr  can  be  obtained  from  Eqs.  (39)  of  our 
previous  report  [1]  in  a  greatly  simplified  form  to  be  written  as  follows; 

tr2p0rV  +  p0R  -  yopaU  +  X*  -  0;  (43) 

<r2po  U  4-  pqR'  +  yopoX  —  po(yo^)’  +  (XJf)'  —  0;  (44) 

r2R"  +  2rR'  —  nOi  +  i )R  -  ^nGr2{p'QU  +  p0X) .  (45) 

For  typographical  convenience,  we  have  omitted  the  subscript  n  referring  to  the  order  of  the  harmonic 
because  our  considerations  henceforth  will  apply  to  any  given  value  of  n. 

In  the  preceding  equations,  primes  denote  derivatives  with  respect  to  the  radial  distance  r,  we 
have  also  set 

rX-  rU'+  2U-  n(n  +  \)V  (46) 

to  represent  the  dilatation  of  the  material.  If  we  differentiate  Eq.  (43)  with  respect  to  r,  subtract  from 
it  Eq.  (44)  and  simplify  the  ensuing  result  by  means  of  the  original  Eq.  (43),  we  reach  the  following: 

U  +~-\x-<T2(Y+rV’-  U).  (47) 

l  Po  I 

In  the  limit,  as  <r  approaches  zero  while  both  U,  V  remain  finite,  we  either  reach  the  adiabatic  condition 

yo  +  "  0,  (48) 

Po 

also  known  as  the  Adams- Williamson  condition,  or  we  must  have  X  —  0.  This  means  that  the  dilata¬ 
tion  must  vanish. 

Pekeris  and  Accad  have  discussed  the  constitution  of  the  liquid  core  by  means  of  a  nondimen- 
sional  stratification  function  /3  (r)  defined  as: 

y°  +  -  yo/90-)-  (49) 

(See  Ref.  5).  This  in  turn,  is  related  to  the  Brunt- Vaisala  frequency  WJ,  according  to 

fi - Kr  N1.  (50) 

Portf 

Excluding  the  existence  of  neutral  stratification  throughout  the  whole  liquid  core,  i.e.,  /3  =  0  every¬ 
where,  we  realize  that,  as  <r  approaches  zero  X  must  also  tend  to  zero.  Soiving  then  Eqs.  (43),  (46) 
and  (45)  under  these  specific  conditions,  we  find  the  particular  solution 
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w  -  R*t y0;  V - r^TTT  lr(u*r  +  2U*] ' 

n(n+l) 

r2(R*)"  +  2 riR*)1  -  [n(n+l)  +  4w Gpor2/y0]R*  -  0. 

We  must  choose  that  expression  for  R*  that  remains  finite  at  the  origin,  r  —  0;  this  will  give  rise  to 
only  one  arbitrary  parameter.  Furthermore,  if  we  allow  that  the  variable  V*  be  discontinuous  at  the 
core-mantle  interface,  we  will  have  available  only  two  parameters  for  the  purpose  of  satisfying  the  three 
boundary  conditions  that  are  valid  at  the  loading  surface,  r  -  a. 

One  plausible  physical  condition  that  can  justify  the  use  of  an  extra  free  parameter  is  the  existence 
of  an  infinitesimal  boundary  layer  at  the  top  of  the  liquid  core,  because  within  such  layer  the  dilatation 
X  can  switch  from  zero  to  a  finite  value  X*.  This  arbitrary  jump  in  value  can  play  the  role  of  the  addi¬ 
tional  free  parameter. 

We  shall  prove  in  a  rigorous  way  that  such  a  condition  really  exists  at  the  core-mantle  boundary; 
we  shall  accomplish  this  by  considering  the  asymptotic  expansion  of  the  equations  of  motion  with 
respect  to  the  inverse  frequency  of  oscillation.  We  shall  limit  our  discussion  to  a  simple  physical  model 
consisting  of  a  uniform  liquid  core  surrounded  by  a  uniform  solid  mantle  so  that  the  analytical  develop¬ 
ments  can  be  simplified  without  detracting  any  essential  characteristic  feature  from  the  phenomenon. 

Let  us  write  the  fundamental  equations  of  motion  for  the  case  of  a  uniform  liquid  core,  that  is 
when  po  and  A  are  constant  and  p.  —  0.  It  is  convenient  to  introduce  the  constant, 

.  4  _ 

A  -  y  IT  Op0, 

which  is  related  to  the  gravitational  acceleration  according  to 

yo(r)  -  Ar, 

and  use  it  to  define  a  new  nondimensional  constant, 

a  —  <r2/A, 

and  a  new  variable  Q  —  R/A,  which  has  the  dimension  of  a  squared  length.  In  this  notation,  Eqs.  (43) 
and  (45)  become: 

~  X-  rW-  a  V)-  Q,  (51) 

PqA 

r2Q"+  2rQ’~  n(n  +  \)Q  -  3r2X.  (52) 

Equations  (46)  and  (47)  can  be  combined  into  a  three-term  relation: 

rX-a(V+  rV-  U)  -  rU’  +  2U-  n(n  +  1)K  (53) 

We  must  solve  Eqs.  (51),  (52),  and  (53)  for  U ,  V,  and  Q.  Let  us  begin  with  Eq.  (53)  and  equate  its 
second  and  third  member;  after  some  algebraic  manipulations,  we  can  write 

r(U  -  aV)'+  (2  +  «)(l/-aF)  -  [n(n  +  2)  -  a(a  +  1)1  K 

It  is  easy  to  realize  that  rl+a  is  an  integrating  factor  for  the  above  equation;  this  will  ultimately  lead  to 
its  solution  in  the  form 

U-  aV-  [n(n  +  1)  -  a(a  +  l)]/(r),  (54) 


where 


r2+°l{r)-  fa  r'+aV(r)dr. 


(55) 
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We  next  express  the  rX  appearing  in  the  right-hand  side  of  Eq.  (52)  by  means  of  the  second  member  of 
Eq.  (53),  and  eliminate  the  U  through  Eq.  (54).  The  right-hand  side  of  Eq.  (52)  can  then  be  written  as 

3ar{rV  +  (1  —  a)  V  —  [n(n  +  1)  —  a(a  +  1)]/). 

The  preceding  expression  is  instrumental  in  establishing  the  fact  that 

3  arl(r) 

is  a  solution  to  Eq.  (52).  This  can  be  verified  by  the  use  of  the  following  relations, 

(r/)'=  V  -  (1  +<*)/; 

r(rl)"  -  rV'  -  (1  +  a)  V  +  (1  +  a)(2  +  a)l\  (56) 

which  are  obtainable  by  differentiating  Eq.  (55)  defining  the  integral  Hr). 


The  general  solution  of  Eq.  (52)  can  be  written  as 

0  -  Br"  +  3 arl(r).  (57) 

where  B  is  the  only  arbitrary  constant  appearing  in  our  formulation  because  0  must  remain  finite  at 
r  —  0.  We  shall  use  Eqs.  (54)  and  (57)  to  eliminate  the  integral  Hr)  between  them  and  express  the 
dilatation  X ,  as  given  by  Eq.  (51),  in  terms  of  0  alone.  This  can  be  used  to  rewrite  Eq.  (52)  as  an 
equation  containing  only  the  unknown  function  0.  If  we  adopt  the  nondimensional  independent  vari¬ 
able 

s  -  r/a. 


the  equation  governing  the  variation  of  Q  becomes 

s 2  -  Q  +  2s  —  n(n  +  1)0  +  — -  [a(a  +  4)  —  n(n  +  I)Js20  —  Ls"+2.  (58) 

dsl  ds  Ka 


Here  L  denotes  the  expression 

-  [n(fl  +  1)  —  «(a  +  \)]Ban, 

Ka 

which  depends  on  the  arbitrary  constant  B.  Equation  (58)  can  be  shown  to  have  the  particular  solution 
Ds ",  where 

n  fl(w  +  1)  —  g(q  +  0  an 
u~  n(n  +  1)  —  a(a  +4) 

Its  general  solution  is  of  the  form 


0  -  Ds"  +  r, 

where  T  is  the  general  solution  of  the  equation  obtainable  from  Eq.  (58)  by  deleting  its  right-hand  side. 
The  latter  equation  can  be  rearranged  into 

■^r(sT)-  IC(s)+p2](sD-0,  (59) 

ds 1 


where 


C(s) 


■n(n-t  0  -  («  +  4)  £2^1, 

A 


and 


u2  mm 


mini  + 1, .  fsd^ 


n(n  +  1)  >  0. 


(60) 
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When  the  frequency  a-,  and  consequently  the  parameter  a,  decrease  to  zero  and  for  s  0,  the  values 
of  the  function  C(s)  will  become  and  remain  negligible  as  compared  with  the  constant  v2.  We 
therefore  conclude  that  for  0  <  s  <  s0,  the  asymptotic  solution  of  Eq.  (59)  valid  for  decreasing  values 
of  <r  is  of  the  form 

r-  -  exp  Ms  -  So) ]  +  0(a) 
s 

where  F  stands  for  an  arbitrary  constant  and  so  —  b/a,  b  being  the  radius  of  the  liquid  core. 

The  asymptotic  expansion  of  Q  accordingly  will  be 

Q  —  Ds"  +  —  exp  Ms  -  s0)l  +  0(a);  (61) 

s 

and  the  asymptotic  expansions  of  the  other  pertinent  variables  U,  K  and  X  can  be  evaluated  therefrom. 


In  fact,  if  we  express  the  first  of  Eqs.  (56)  in  terms  of  the  nondimensional  variable  s,  we  get 

P(s)  -  s  ~  +  (2  +  a)I(s). 
as 

The  expression  for  /(s)  is  obtainable  by  equating  Eq.  (61)  to  Eq.  (57)  since  they  both  represent  the 
Q- variable: 

Q  -  Dsn  +  -j  exp  Ms  -  s0)]  -  Ba"s "  +  3asa/(s). 

Proceeding  along  these  lines  and  neglecting  the  positive  powers  of  the  small  parameter  a  (or  a),  and 
since  we  are  considering  asymptotic  expansions,  we  find,  after  elementary  operations,  that 

F(s)  -  —(as)"-1  +  exp  Ms  -  s0)]. 
n  3  aas 

Similarly,  by  using  Eqs.  (51)  and  (54),  we  shall  reach  the  asymptotic  representations  of  U  and  of  the 
normal  stress  \X\  the  results  are: 

Uis)  -  Bias )"-1  +  exp  Ms  —  s0)l; 

3aa  si 

Xis)  -  -5-  —  exp  Ms  -  s0)J. 

3  <r  s 


If  we  introduce  the  nondimensional  arbitrary  parameters 

E\  -  Ban~2,  F]  -  F/3aa2, 
we  can  write  the  previous  results  in  dimensionless  form  as  follows: 


-  P(s) 

a 

-  Uis) 
a 

Ris) 

ayoia) 

Xis) 


—  E\  s"-1  +  —  i iF]  exp  Ms  -  s0)]; 
ft  s 

£,s"-1  +  — -  -  y-  —  F i  exp  Ms  -  s0)]; 
sl 

_  £<jL  -  [i  -  «-|£ +  bfl  exp  Ms  -  s0)l; 

a1  A  a2  (  n  j  s 


F i  exp  Ms  -  s0)]. 


(62) 


These  equations  prove  that  within  the  liquid  core  the  oscillations  induced  by  a  long-period  forcing  func¬ 
tion  vary  as  exp  (1/cr )  and  not  linearly  with  or2  as  one  might  have  originally  surmised  by  inspection  of 
Eqs.  (43)  and  (44). 
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The  monomial  terms  that  appear  in  the  expressions  for  U  and  R  are  regular  functions  in  the  inter¬ 
val  0  <  i  <  s0  and  will  offset  any  discontinuity  that  might  arise  because  of  the  presence  of  the 
exponentials.  The  F  variable  varies  as  v  —  Mo  and  will  tend  to  increase;  we  know,  however,  from 
other  sources  that  F  will  be  discontinuous  at  the  interface. 


On  the  other  hand,  the  variable  X  related  to  both  the  dilatation  and  the  normal  stress  in  the  liquid 
core  depends  on  the  exponential  alone  since  <*v2  is  a  finite  constant;  the  X  will  therefore  exhibit  a 
discontinuous  behavior  typically  induced  by  a  boundary  layer  as  described  in  the  following. 

The  behavior  of  the  exponential  term  depends  on  the  relative  values  of  (s  —  s0)  and  v  Mo. 
For  small  values  of  o  and  for  s  much  less  than  %  s  ~  so  >s  negative,  and  the  exponential  of  a  large 
negative  number  provides  a  negligible  contribution;  we  can  say  then  that,  from  a  practical  point  of 
view,  the  X  variable  vanishes  at  large  distances  from  the  top  layer. 


On  the  other  hand,  consider  the  situation  where  s  -  s0  is  small  enough  to  be  comparable  in  mag¬ 
nitude  with  the  assumed  small  value  of  o ;  the  exponential  then  provides  a  sizeable  contribution.  The 
smaller  the  value  of  o ,  the  thinner  the  width  of  the  topmost  layer  where  X  is  nonvanishing.  In  the 
limit  as  a  approaches  zero,  the  stress  distribution  can  be  assumed  to  behave  as  a  delta-function  having 
a  finite  spike  at  $  —  So- 


For  s 


so 


b/a,  Eq.  (62)  become 

-  U(b)  -  E\  sS~'  +  -(~t  1-  F,; 
a  stf 

X{b )  -  —  F,  -  n(n  +  1)  F,; 

So  Xs0 


R(b) 

ay0(a) 


E\  sfii 


(63) 


whereas  Fean  be  taken  to  be  discontinuous  at  the  interface.  We  have  two  arbitrary  parameters  Eh  F, 
plus  the  value  of  the  discontinuity  suffered  by  Fin  going  from  the  liquid  to  the  solid  layer.  These  are 
then  three  quantities  we  can  use  in  trying  to  satisfy  the  boundary  conditions  at  r  —  a. 


This  approach  was  followed  in  evaluating  the  load  numbers  of  low  order,  e.g.,  of  order  n  =■>  2,  3, 
4,  which  are  rather  well  known  from  direct  physical  measurements.  The  agreement  was  good.  We  are 
therefore  convinced  that  our  analytical  solution  is  an  adequate  representation  of  the  physical  conditions 
existing  at  the  liquid  core-solid  mantle  interface;  and  believe  that  the  same  procedure  should  be  used 
for  the  calculation  of  higher  order  load  numbers  which  are  not  known  from  ground  data  but  which  can, 
nevertheless,  play  a  significant  role  in  the  final  determination  of  the  spheroidal  deformations. 
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